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Abstract 

From the study of the multiple muon events in deep underground detectors, it is 
possible to extract information about the spectrum and composition of the primary 
cosmic rays. In this work the number of TeV muons produced by a primary cosmic 
ray of a given energy and zenith angle is computed using analytic and montecarlo 
methods, for a family of simplified models as description of the properties of hadronic 
interactions. The effects of the uncertainties in our knowledge of the hadronic cross 
sections in the calculation of TeV muons are discussed. 
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1 Introduction 



Deep underground detectors measure events with several nearly parallel muons fl], 0, 0] , 
separated by few meters 0]. The muons are produced by the decay of charged pions and 
kaons in the hadronic showers induced by primary cosmic rays. From the study of these 
events it is possible to obtain information about the spectrum and composition of cosmic 
rays in the energy region (10^^ ^ Eq ^ 10^^ eV). 

The shower produced by a nucleus of total energy Eq and mass number A can be de- 
scribed with reasonable accuracy as the superposition of A proton showers each of energy 
Eq/A. Therefore if a proton of energy Eq produces (above a threshold energy -Emin, and 
with a zenith angle dependence oc (cos6')~^) an average number of muons {N^{Eo))p, a 
nucleus of mass number A and the same total energy will produce with good approxi- 
mation an average number of muons: {N^{Eo))a — A{N^{Eq/ A))p. The average 
muon multiplicity of a proton shower grows more slowly that the primary energy, approx- 
imately as E^''^^; this can be understood observing that in the approximation of Feynamn 
scaling the number of mesons above energy -Emin in a shower, neglecting threshold effects, 
grows linearly with energy, but with increasing energy the mesons are produced deeper 
in the atmosphere, where their decay is more rare. Therefore the average number of high 
energy muons in a shower of energy Eq depends on the mass number of the primary nu- 
cleus approximately as oc A^'"^^, and (for the same energy spectrum) a primary comic ray 
flux rich in heavy nuclei will produce more high multiplicity events in deep undergroud 
muon detectors than a flux with a 'light' composition. This allows us in principle to ob- 
tain information on the primary cosmic ray spectrum and composition from underground 
measurements of multiple muon events. The uncertainty on the properties of the hadronic 
interactions is the main source of systematic errors in the development of this program. 

In this work we will discuss some methods to estimate the size of this uncertainty. 
Elbert |^ has suggested that the average number of muons above a minimum energy £"111111 
produced in the shower of a primary proton of energy Eq has the scaling form: 

(iV,(£.ni„, Eo)) = ^ ^ , G (^) , (1) 



COS 9 V Eq 

Gaisser and Stanev ^ and Forti et al. have also fitted their montecarlo results with 
the form (0). In section 2 we will derive formally this result, discussing under which 
conditions it is valid, we will also show how it is possible to compute with analytic 
methods the function G{x) from a knowledge of the inclusive single-particle differential 
cross sections. 

In section 3 we will discuss a family of 'toy models' for the hadronic cross sections that 
are a generalization of the algorithms introduced by Hillas [ffUl. These models are fully 



described by a small set of simple montecarlo algorithms, and are therefore constructed as 
montecarlo instruments, at the same time the inclusive differential cross sections generated 
by the algorithms have simple closed form analytic representations, and are also suitable 
for analytic studies. 

In section 4 we will discuss some explicit calculations of high energy muon produc- 
tion in proton showers, comparing analytic and montecarlo methods. The montecarlo 



2 



technique is necessary for a correct treatment of fluctuations in the shower development. 

In section 5 we will compare two models of hadronic interactions that are constructed 
to produce the same inclusive muon flux, with one model having a larger and softer vr 
multiplicity, and discuss how the multiplicity distribution of TeV muons changes. 



2 Analytic Method 

In this section we will discuss how under two approximations: (i) energy independent 
interaction lengths, (ii) validity of Feynman scaling in the fragmentation region of the 
hadronic cross sections, it is possible to solve analytically the shower development equa- 
tions and calculate the inclusive high energy muon spectrum produced by a primary proton 
of energy Eq. The solution has the scaling form (||) suggested by Elbert and used in 
P, ^, (there is also a differential form of this scaling law given in equation (P)). The 
function G{x) depends on the interactions lengths of nucleons, and mesons, and on the 
inclusive differential cross sections daa^b/dx for production of particle a in the interaction 
of particle b with x = Eb/Ea. There is also a dependence on the density distribution of 
the medium where the showers develop. For an exponential distribution: p{h) oc e^'*/''", 
as is the case for cosmic rays reaching the earth atmosphere and zenith angles not too 
close to the horizontal direction, this results in G being proportional to the scale height 
ho- 

The analytic calculation of the functions G can be performed with a straightforward 



generalization of the methods used for the calculation of electromagnetic showers [O. 
We will begin to discuss the inclusive muon flux produced by a power law flux of primary 
particles, and then discuss the case of a monochromatic beam of particles. 

2.1 Power law initial proton spectrum 

The calculation of the inclusive muon flux produced by a power law primary flux, is a well 
known problem (see for example the textbook |^), and is directly applicable to the real 
flux of primary particle reaching the earth. The all nucleon cosmic rays spectrum is in 
fact well represented by a power law: 0o(-Eo) — K Eq" with index a = 2.7, and a constant 
K ~ 1.85 (with Eq in GeV and in (cm^ s sr GeV)~^). Using the approximations of 
constant interactions lengths and Feynman scaling, the muon spectrum above an energy 
of ~ 20 GeV has the approximate form [H] : 



L,r(«) Ecos9 



-1 



LkM Ecos9 
1 H 



-1 




We will be mostly interested in the high energy limit {E ^ e^,ej^): 

^ e.HAa) +eM) ^ _ jM. K E- (3) 

E cos 9 E cos 9 
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The critical energy = rriT^hQ/cTT^ ~ 115 GeV can be interpreted approximately as the 
energy at wich the average charged pion produced in a shower has a probability of 1/2 to 
decay, Ek — 850 GeV has the same meaning for charged kaons, {rriK) and r^r {tk) are 
the mass and lifetime of a charged pion (kaon), = 6.34 km is the scale-height of the 
stratosphere. The 'constants' L.^, L^, H^^ and H^, depend on the esponent of the primary 
flux a, and on the properties of hadronic interactions. Explicit solutions for Lj^{a) and 
HJa) are 



HJa] 



1 



- A 



In 



N 



a(l 

A, 
A 



N 



1 - r: 



a+1 



(« + l)(l-r^) 



In these equations: = {m^/mT,Y] the quantities Zab{ci) are defined as: 



Zabioi) 



dx X 



a-l 



' dn{x) 
dx 



(4) 
(5) 

(6) 



where dua^b/dx is the spectrum of particle b in the interaction of particle a with a target 
air nucleus and x = Eb/Ea, A^r = AAr/(l — Ztvat), A^^ = A^/(l — Z^rTr) where A^v and A^^ are 
the interactions lengths in air for nucleons and pions. The constants Lk{o) and Hk{c() 
have similar expressions, the inclusion of the processes (vr^ — > K^) and {Kl — > K^) 
introduces some small complications |T^. The flux of only positive or negative muons can 
also be calculated with expressions similar to (H). 



2.2 Monochromatic initial proton spectrum 

The flux of muon produced by a monochromatic beam of protons takes a more complicated 
form. We can define the function f{E;EQ), the inclusive differential spectrum of /i's 
produced by a primary proton of energy Eq, and g^E^i^i; Eq) the integral spectrum (or 
average number of muons) above energy -Emin- The functions / and g are related by: 



g{E^^- Eo) = {N^iE^,^, E,)) = / dE f{E- E^) (7) 
If E (-Emin) ^ ^TTi^Ki equation @ implies: 



roo 

/ dEoE^" f{E,Eo) 



= Ua) = HJa) + Ek HKia) 



This condition can be satisfied only if the functions / and g have the scaling form : 

mEo) = -^^^,F(^) (9) 



cos^E2 y^^ 
1 
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Equation (^) defines also the Mellin transforms Mp and Mq of the function F and G: 



M^(s) = e^(s + 2) (11) 

Ma{s) = ^^i^^ (12) 
The Mellin transform of a functions A{x) defined in the interval < x < 1 is: 

Ma{s) = C dx x'A{x) (13) 

From equation (|^) it is possible to deduce the following relation between the functions F 
and G : 



F{x) = -x^4- 
dx 



G{x) 



X 



(14) 



that implies in general Mg{s) = Mf{s)/{s + 4). 

For continuous A{x) the transform Ma{s) is well defined in the complex field for all 
values of the imaginary part of its argument and for the real part in the interval 

Si < 3?[s] < 52- The relation ( [T^ ) can be inverted [T^, [TT| with: 

Aix) = — I ds Ma(s)x~'^'^'^ (15) 
2m Jc 

where the integral is taken in the complex field along an arbitrary path C that starts 
from Q{s) = — oo and ends at Q'(s) = +oo, and is entirely contained in the domain of 
definition of Ma{s). 

To summarize: the function e^(a), introduced in (^, for a real and positive determines 
the high energy inclusive muon flux produced by a power law primary flux of index a: 
4>^{E)/(I)q{E) = e^{a)/{E cos 9). The function e^(a) can be calculated from the inclusive 
hadronic cross section according to equations and (^, and the definition is valid also 
for a in an infinite region of the complex field. According to equations (|ll]) and (0), 
e^(a) gives also the Mellin transform of the functions F{x) and G{x) that determine the 
differential and integral inclusive muon spectrum for a monochromatic beam of protons. 
We can therefore compute the functions G{x) and F{x) in two steps: (i) compute the 
function e^(a;), (ii) use the inversion formula (|T3|). 



3 Models for the differential cross sections 

In this section we will discuss an explicit family of models for the hadronic cross sec- 
tions, that is a simple generalization of the montecarlo algorithms proposed by Hillas 
[ Pmi , to model multiple particle production in high energy hadronic interactions. These 
algorithms, in spite of their simplicity, give results in good agreement with data. The 
inclusive differential cross sections generated by the montecarlo algorithms have also ex- 
plicit and exact analytic representations, this allows to cross check results obtained with 
analytic and montecarlo methods. 
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To explore the sensitivity of muon production to the characteristics of the cross sec- 
tions, we have generahsed the set of algorithms originally proposed by Hillas, introducing 
some parameters that can be changed continuously. In chapter 17 of reference 0, Gaisser 
describes in some detail the Hillas algorithms, and discuss also possible generalizations. In 
[|T0[] Hillas considered only nucleons and pions, neglecting kaons and other particles. It is 
rather straightforward in the framework of the model to include kaons or treat separately 
protons and neutrons, but for the sake of simplicity we will also consider only nucleons 
and pions. The separate treatment of protons and neutrons is of little importance if we 
sum over the charge of the muons; only ~ 25% of TeV muons are produced by kaons, and 
therefore, for our illustrative discussion, kaons are not of crucial importance. The main 
focus of this work is not to obtain an absolute calculation but to discuss the sensitivity 
of the results to variations of the 'input'. The algorithms that we are considering are 
exactly scaling in the variable x = E/Eq where Eq is the incident particle energy. Scal- 
ing violations can be be introduced allowing the parameters of the model to vary with 
energy, we will however discuss here only scaling models, and concentrate on the study 
of the effects of distortions of the spectra of nucleons and pions. We will not discuss the 
transverse momentum of the particles produced in the shower, and we will not comment 
on the separation between muons. 

3.1 Proton interactions 

The set of algorithms to generate the interaction of a proton of initial energyy Eq are: 

1. Generate x (0 < a; < 1) with linear distribution of average (x) = 1 — Kp. Construct 
a leading nucleon of energy Eqx. 

2. The remaining energy Eo{l — x) is divided into two parts A and B with a flat 
distribution. Then each piece of energy is split again into two parts with a flat 
distribution {A ^ Ai + A2 and B Bi + B2). 

3. Each of the four pieces of energy has now a probability P* of splitting once again 
into two parts, always with flat distribution. 

4. At this point we have N energy pieces (4 < < 8), with a binomial probability 
distribution p{N) = Pbinomiai(A^ — 4; 4, P*). Each piece of energy now enters a 
recursive process. 

5. The energy of the piece is splitted with flat distribution into two parts. A part is 
assigned as the energy of a pion (with equal probability for the three charges), the 
remaining energy is again divided into two parts: one is assigned to a new pion, 
for the other the process is iterated. The splitting is stopped when the remaining 
energy is less than a pion mass or of some preassigned threshold value. 

In the algorithms we have introduced two parameters: the inelasticity Kp (1/2 < Kp < 
2/3) that gives the fraction of the projectile particle energy that is not transfered to the 
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leading nucleoli, and P* (0 < P* < 1) that determines the shape of the pion spectrum, 



when (x — > 1) the pion spectrum has the behaviour oc (1 



■ X 



i3 ; 



if P* = and oc (l-x)^ if 



P* = 1. A 'model' in this framework is then represented by the pair of numbers {Kp, P*}. 
The original algorithms of Hillas corresponds to the model {|, 0}. 

The montecarlo algorithms described above give inclusive differential cross sections 
that have explicit analytic expressions: 



(16) 



dx 



1-6 K, 



3 

96 / ^ 1 



1 ^ , Inx 2 
1 + lnx- ^ — 

X 2 



'l + P*) + 2P*^i^^| + 
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(17) 



The differential cross section for 7r° production is simply 1/2 of (p!7|). The plateau of 
the rapidity distribution for the charged pions has a height that increases linearly with 
(1 + P*): 



Pp—f-K^ 



-[1 + P*] 



1^ 



The momenta of the inclusive spectra Zpxia) = (s" )px can be calculated explicitely: 



1-P* 



1 + Q K 



--12(k„-^ , 
a V 27 a + 1' 



a - 2 



a 



1 / _ / 12 _ 6 



(19) 



(20) 



For a = 2, Zpp = 1 — Kp and ZpT^± = ^Kp independently from P* and the shape of the 
spectra. 



3.2 Charged pion interactions 

To generate a vr^ interaction the algorithms are the following 

1. The total energy is divided into two parts A and B with a flat distribution. The 
first part A is assigned to a pion with probability Po- This 'diffractive' pion is 
charged with a probabihty of 0.87, the probability of being of same (opposite) sign 
with respect to the projectile is 0.80 (0.07). 

2. The energy of the piece B is divided into two parts Bi and B2, then each piece is 
again divided into two. We now have 4 energy pieces. Bn, B12, P21 and P22- 

3. Piece A if it was not already assigned to the diffractive pion, is treated in the same 
way as piece B in the previous and following steps. 
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4. The pieces Bu and B21 have now each a probabihty Pa of being assigned to a pion 
(with equal probabihty for the three charges). 



5. Each of the remaining energy pieces (-B12, -B22 and Bu, B21 if not aheady assigned 
to pions) are sphtted into two parts with flat distribution with a probability Pb 

6. Part B is now divided into (2 < < 8) energy pieces. Each of these pieces is 
fragmented into particles, with a recursive process (as in proton interactions). The 
energy of one piece is sphtted with flat distribution into two parts, a part is assigned 
as the energy of a pion, for the remaining energy the process is iterated until we are 
left with a remaining energy less than a pion mass or a preassigned threshold value. 

The pion interaction model depends on three parameters {Pd, Pa-, Pb}- The original algo- 
rithm proposed by Hillas [1^ corresponds to the choice {|, |, 0}. The physical meaning of 
the parameters is easy to understand intuitively. Increasing and Pa the pion spectrum 
is hardened and the multiplicity decreased. Increasing Pb has the opposite effect: the 
spectrum is softened, and the multiplicity increases. 

The inclusive pion spectrum implied by this set of algorithms can be calculated: 

0.87 Pd + - (2- Pd)Pa^^+ (21) 



dx 3 

+ i(2-P^)(2-PD) 








The spectrum of 7r° is obtained changing the numerical coefficients of the three terms in 
(|1]): 0.87 0.13, 4/3 2/3, and 4/3 2/3. The momenta of the energy distribution 
are: 



= + 1(2 - Pn)PA^ + -A2 - PA)i2 - Pb) ^ 



a 3 3 a^{a — 1] 

The rapidity density is: 



a 



(22) 



p^^ = ^(2 - P^)(2 - Pb)(1 + Pb) (23) 



4 Explicit Calculations 

In this section we will discuss some calculations of TeV muons, using the algorithms 
discussed in the previous section to describe hadronic interactions. We will always assume 
that the interactions length are constant and that Feynman scaling is valid. In our 
framework therefore an 'hadronic interaction model' is fully described as a 'vector' of eight 
numbers: M = {Xp, A^; Kp, P*; Pb, Pa, Pb}, the first two quantities are the interactions 
lengths of nucleons and pions, the next two numbers describe protons interactions, the 
last three numbers describe pion interactions. 
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The first quantity that is interesting to calculate is the inclusive high energy muon 
flux. This flux is determined by the quantity e^(a) for a = 2.7. In table 1 we show the 
value of e^(2.7) calculated for different 'hadronic interaction models'. 



Table 1. e^{a = 2.7). The primary flux has form: (polE) = K E the muon flux 
{E > 100 GeV) has form: = e^/{E cosO) x ^E). 



Xp (g cm-2) 


(g cm-2) 


Kp 


p* 


Pd 


Pa 


Pb 


(GeV) 


86.0 


111.8 


1/2 





1/2 


1/2 





8.57 


77.4 


111.8 


1/2 





1/2 


1/2 





9.02 


86.0 


100.6 


1/2 





1/2 


1/2 





8.14 


77.4 


100.6 


1/2 





1/2 


1/2 





8.57 


86.0 


111.8 


1/3 





1/2 


1/2 





5.42 


86.0 


111.8 


2/3 





1/2 


1/2 





11.03 


86.0 


111.8 


1/2 


1 


1/2 


1/2 





6.35 


86.0 


111.8 


1/2 








1/2 





7.96 


86.0 


111.8 


1/2 





1 


1/2 





9.34 


86.0 


111.8 


1/2 





1/2 








8.45 


86.0 


111.8 


1/2 





1/2 


1 





8.70 


86.0 


111.8 


1/2 





1/2 


1/2 


1 


8.44 



The flrst line of table 1 corresponds to the differential cross sections proposed originally 
by Hillas [T^, and is also our choice for a 'reference model', the other lines show the 
sensitivity to 10% changements in the interactions length, and to the maximum allowed 
variations of each of the parameters that we have constructed to describe the shape of the 
differential cross section. We can observe that the inclusive muon production is especially 
sensitive to the properties of the proton interactions. Decreasing the inelasticity Kp or 
softening the pion spectrum (increasing P*) depresses the muon flux. Modiflcations of 
the pion differential cross section have a smaller effect on the inclusive muon flux. In fact, 
because of the steepness of the primary flux, most of the muons are produced in the decay 
of mesons produced in the flrst interaction of a primary cosmic ray. 

The function G{x) that gives the average number of high energy muons produced by 
a primary particle is easily calculated from (0) using the inversion formula (|15[). In the 
models we are considering ep{a) can be obtained from (||), (|^), (|l^), (|20D , (^21) , as a simple 
combination of elementary functions, and the Mellin transform can be inverted with an 
easy numerical integration. 

The function G{x) calculated as discussed above for the 'reference model': [ Ap = 86 
and A^ = 111.8 g cm~^, {Kp,P*] Po,Pa,Pb} = {|;0; |,|,0}], is shown in flgure 1. 
In the same flgure we also show for comparison the curves that Gaisser and Stanev 
and Forti et al. have used as flt to their montecarlo calculations. The agreement is 
surprisingly good considering the extreme simplicity of the model we are discussing. It 
should also be noted that the results of [|, refer to muons not above a flxed threshold 
energy, but at a flxed depth h. For the comparison we have used the approximation 
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Emin{h) 0.53 (e^-^^ — 1) {h in km.w.e., E in TeV), fluctuations in the muon energy loss 
should however be taken into account for a more detailed comparison. 

In figures 2,3,4,5 we show how the function G{x) is modified because of changements 
in the shape and normalization of the inclusive differential cross sections. We have re- 
calculated G{x) using different values of the parameters. We will consider the 'reference 
model', and change one parameter at the time. 

In figure 2a we show the {p — > vr^) spectrum obtained with different values of the 
inelasticity Kp = |, |, | and P* = 0. In figure 2b, we show the function G{x) as 
a function of 1/x = E^/E^i^ for these three values of Kp. As expected the number of 
muons increases for larger Kp, because more energy is transfered to pions. With increasing 
Eq/E^[^ however the difference becomes smaller, the curves join, and in fact on close 
inspection cross each other. This can be explained observing that for large i?o/-Emin it 
is possible to obtain more pions above the threshold energy -Emin giving more energy 
to the leading nucleon in the interaction. The number of pions in the first interaction 
of the shower decreases, but additional pions are created in the second interaction, our 
calculation takes into account the fact that these additional pions are produced deeper in 
the atmosphere and have a smaller probability to decay. 

In figure 3a we show the {p n'^) spectrum for P* = 1,0 {Kp = |). As discussed 
before P* = 1 corresponds to a softer spectrum, with a higher multiplicity. In figure 3b 
we show the function G{x) as a function of 1/x = Eq/E^^^ calculated for the two pion 
spectra. For small i^oZ-E-min the curve corresponding to P* = (harder vr spectrum) is 
higher; the two curves cross each other at i?o/-Emin — 30. This is qualitatively easy to 
understand: for small £'o/-Emin the harder spectrum produces more pions above threshold, 
when Eq/ E^in grows, the effect of the larger pion multiplicity becomes dominant. 

In figure 4a and 5a we show the (vr^ — > vr^) spectrum for different values of the 
parameters : Po = 0, |, 1, and Pg = 0, 1. In figure 4b and 5b we show the curve G{x) 
calculated with the different spectra. Some remarks are : (i) very large deformation of 
the (vr^ TT^) spectrum result in small variations of {N^^E^i^, Eq)); (ii) the effect is 
very small for Eq/E^i^ < 10 when most muons come from the decay of first generation 
mesons; (iii) the largest effects comes from different treatments of vr diffraction (fig 4b); 
(iv) a softening of the vr — > vr spectrum from dn/dx oc (1 — x)^ [Pb = 0) to oc (1 — xY 
{Pb = 1) depresses muon production for Eo/Emm ~ 10^, then enhances it. The effect (see 
fig 5b) is however small. 

To summarize the information of the effect of G{x) of the changement of the parameters 
of the model, it can be useful to discuss the quantity .^p(-Emin/ Eq), the logarithm derivative 
of the average muon multiplicity as a function of the parameter P taken from the 'starting 
model' that we take as the original Hillas model: {Kp, P*; Pd, Pa, Pb} = {|)0, |,|,0} 
with Xp = 86 and = 111.8 g cm~^. 



The meaning of C,p is that if the parameter P is changed by (for example) 10%, the 




1 



(24) 
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resulting percentual effect on {N^^E^ii^, Eq)) is (^p x 10)%. A positive (negative) C,p 
indicates that an increase of parameter P will produce a depression (enhancement) of 
(A^"^). The logarithimic derivatives C,p{x) are shown in table 2. 



Table 2. Logarithmic derivatives ^p(x). (x = Eram/Eo). 



Parameter / x 


10-1 


10-2 


10-3 


10-^ 


10-5 


Xp 


-0.530 


-0.497 


-0.411 


-0.333 


-0.276 




0.512 


0.483 


0.405 


0.333 


0.280 




0.985 


0.230 


0.082 


0.033 


-0.002 


(1 -P*) 


-0.282 


0.133 


0.192 


0.173 


0.155 


Pd 


0.082 


0.107 


0.077 


0.041 


0.012 


Pa 


0.012 


0.033 


0.004 


-0.043 


-0.088 


(I-Pb) 


-0.009 


-0.053 


-0.025 


0.054 


0.131 



The first two lines of table 2, show the dependence on the hadronic interaction lengths. 
An increase in Xp decrease the number of muons because with a larger proton interaction 
lengths the shower develops deeper in the atmosphere where the density is higher and 
muons decay is more difficult. When i?o/-E'min grows the effect of a changement of Xp 
decreases in importance because a growing fraction the muons is produced in a cascade 
of type p ^ 71 ^ n —>■ fi. 

An increase of Att increases the number of TeV muons because if the interactions length 
is longer the pions have more time to decay. The effect becomes smaller with increasing 
Eo/E^iri, because a larger A^r will also produce a deeper shower, and the pions produced 
in a cascade of type p ^ vr — > vr are created at lower altitude and decay more rarely. Note 
that the two hadronic interaction lengths enter in the expressions for muon production 
only in the combination: A^r/Ap ~ X-^/Xp and therefore ^(A^r) — —^{Xp). 

The other lines in table 2 describe the dependence on the differential cross sections, 
the same comments developed for the discussion of figures 2-5 apply, one may notice the 
changements of sign of ^{Kp) (at very large Eo/E^i^) and of C,{P*). It is encouraging 
to see that |^p| is always less than 1, showing that the sensitivity to distortions of the 
spectral shape is only moderate. Note how G{x) is relatively insensitive to the properties 
of pion interactions. 

4.1 Fluctuations 

The analytic method that we have described allows us to compute the average number 
of high energy muons in a shower, but does not take into account fluctuations in the 
shower development. In order to study the importance of fluctuations and also to check 
the analytic calculation we have prepared a montecarlo implementation of the family of 
interaction models discussed in section 3. The straightforward montecarlo method is based 
of the shower code developed in Bartol [|r3| (see also |^) and is based on the following 
steps: (i) a primary cosmic ray is propagated in the atmosphere until it interacts; (ii) 
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a set of secondary particles (pions and nucleons) are produced at the interaction point 
according to the algorithms described in section 3; (iii) each one of the secondary particles 
is propagated until it interacts or decay; (iv) at each decay or interaction vertex the 
incident particle is destroyed and a set of new particles is produced conserving energy 
and momentum; (v) the procedure is iterated until all particles are below a preassigned 
minimum energy. All produced muons are recorded. 

In figure 6 we show the inclusive muon spectrum above 1 TeV produced by a vertical 
primary proton of energy 10, 10^, 10^ and 10^ TeV, obtained with the montecarlo method 
and with the analytic formula using our 'reference model'. The results obtained with 
the two methods are in excellent agreement between each other. The montecarlo method 
allows to study not only (N^j) but also the probability distribution P{N^) of having exactly 
Nfj_ muons in a shower. The distributions P{Nf^) calculated with the montecarlo method 
for Eq = 10 and 10^ TeV {Ef'^ = 1 TeV, ^ = 0° ) are shown in figure 7a and 7b. In 
the same figures we compare the montecarlo results with a poissonian distribution of the 
same average and a negative binomial distribution of same average and dispersion. The 
distribution P{N^) is broader than a poissonian, and the difference becomes more marked 
with increasing energy, the negative-binomial being a good fit. These results have been 
found previoulsy by Forti and collaborators 0; we would like to stress that the non- 
poissonian fluctuations are not connected to violations of Feynman or KNO scaling, and 
are present also in exactly scaling models as those we are discussing. 

5 Interaction Model and Multiplicity Distribution 

As an illustration of the importance of the interaction model for the calculation of the 
multiplicity distribution of underground muons, we have calculated the fluxes of TeV 
muons using a 'realistic' proton flux of energy spectrum: (f)Q{E) = 1.85 i?"^ '', steepening 
to oc E~^ ioT E > 3 X 10^ GeV {E in GeV and (f) in (cm^ s sr)~^) ), and two different 
models for the hadronci interactions. In both models Xp = 86 and A^r = 111.8 g cm~^, 
and the charged pion interactions: are described by the algorithms originally proposed by 
Hillas with parameters {Pd, Pa, Pb} = {^,^, 0}. The two models differ in the treatment 
of proton interactions. The first model is nearly identical to the 'reference model' with 
parameters: {Kp, P*} = {0.475, 0}; the second model has a larger inelasticity but a softer 
spectrum: {Kp, P*} = {§, 1}- The G functions obtained with the two models is shown in 
figure 8, both models have e^(2.7) = 4.7 Mg'(0.7) ~ 8.16 GeV, and for energies E ^ e^r 
produce essentially identical 'inclusive' muon fluxes. The larger multiplicity of model-2 
exactly compensates its softer spectrum. We may however expect that using model-2 the 
probability of having several muons in the same shower is larger. 

To investigate quantitatively this possibility we have generated approximately 3 million 
vertical showers with Eq > 1 TeV for each of the two models. In order to increase 
the statistics of events with high muon multiplicity we have sampled the energy of the 
showers from a distribution oc i?o"^'^^5 weighting each event with E^''^^ / (pQ^E). In figure 
9 we show the obtained inclusive muon fluxes that are essentially identical in the two 
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models, and very well represented by: (f)fj,{E) = e/E (poiE) with e = 8.16 GeV. The muon 
multiplicity > 1 TeV) for the two models is shown in figure 10a, and the ratio of 
the fiuxes obtained the two models for the same multiplicity is shown in figure 10b. The 
inclusive fiuxes (pi + 202 + . . . + ncpn + ■ ■ ■ are equal to better that 1%, but the fiux of 
single muons is 5.5% smaller using model-2, the ratio model-2/model-l becomes 1.22 for 
double muons, grows to 1.46 for triples, to 1.76 for quadruple muons, and then seems to 
remain approximately constant. 

The two models considered are 'physically consistent', both respect conservation laws 
(energy, momentum, baryon number), and both produce the same inclusive muon fiux, 
however the same experimental multiplicity distribution of underground muons if inter- 
preted with model-1 (model-2) would result in a heavier (lighter) composition, because 
the effects of the smaller (larger) frequency of high multiplicity events should be compen- 
sated with a different mass distribution of the primaries. 

In this work we do not attempt a more realistic and complete discussion, that will be 
presented in a future paper. We note that the fairly extreme distortions of the spectra 
that we have tried, produce a difference of about a factor of 2 for the frequency of events 
with multiplicity > 10. The present range of uncertainty in composition ||^ can produce 
larger differences. 



6 Conclusions 

In this work we have formally derived the result that the average number of high energy 
muons (-Emin ~ 1 TeV) produced by a primary cosmic ray proton of energy Eq has the 
scaling form: {N^{E^i^, Eq)) = G{E^i^/ Eq) / E^^^. The function G{x) is calculable ana- 
lytically from a knowledge of the inclusive single-particle differential cross sections. We 
have illustrated how the shape and normalization of G{x) depends on the detailed form 
of these cross sections. 

We do expect detectable deviations from the scaling behaviour. A source of deviation 
is simply the fact that the critical energy for kaon decay Ek is not small with respect 
to 1 TeV. A second source of deviation is due to the fact that the hadronic interactions 
lengths are decreasing with energy. There is also the possibility of observable violations of 
Feynman scaling in the fragmentation region, the measured growth of the central plateau 
should have visible effects for large Eo/E^i^. 

We have discussed a possible generalization of the montecarlo algorithms originally 
developed by Hillas |0, to describe the properties of of hadronic interactions. These 



algorithms because of their remarkable simplicity and flexibility, can be a useful too to 
study in detail the effects of uncertainties in the properties of hadronic interactions in the 
development of showers. They could be very useful for the study of the highest energy 
cosmic rays {E ~ 10^° eV). 

Uncertainties in the modeling of hadronic interactions are the dominant source of 
systematic error for the measurement of cosmic ray composition from data on multiple 
muon events in deep underground detectors. The spectrum of the leading nucleon, and 
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of fast pions produced in nucleoli interactions are of special importance. The details of 
particle production in pion interactions are less important to control. 
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;ure Captions 



Fig. 1. Curve G{x) (equation calculated with the Hillas model. The curve 
is compared with the parametrization of the average number of muon at depth h 
of Gaisser and Stanev and Forti et al. assuming the approximate relation 
Emin = 0.53 (e°-^^ - 1) {E in TeV, h in km.w.e.). 

Fig. 2a. Plot of the inclusive (p — > vr^'') spectrum calculated with P* = and with 
inelasticity Kp = |, |, |. 

Fig. 2b. Curve G{x) calculated with the 'reference model' [{Kp,P*; Pd, Pa, Pb} = 
{|,0, |, |,0}, Xp = 86, = 111.8 g cm~^] and with modified inelasticity Kp = |, 

2 

3' 

Fig. 3a. Plot of the inclusive {p —>■ tt^) spectrum calculated with P* = and 
P* = liKp = i). . 

Fig. 3b. Curve G{x) calculated with the 'reference model' (P* = 1) and with the 
modification P* = 1. 

Fig. 4a. Plot of the inclusive (vr^ vr^) spectrum calculated with the model 
{Pd, Pa, Pb} = {Pd, |, 0} and Pd = 0, i, 1. 

Fig. 4b. Curve G{x) calculated with the 'reference model' (Pd = |) and with 
Pd = 0, 1. 

Fig. 5a. Plot of the inclusive {n^ tt^) spectrum calculated with the model 
{Pd, Pa, Pb} = H |, Pb} and Pb = 0, 1. 

Fig. 5b. Curve G{x) calculated with the 'reference model' (Pb = 0) and with 
Pb = 1. 

Fig. 6. Plot of f{Ep_;Eo) the inclusive differential muon spectrum produced by a 
vertical primary proton for Eq = 10, 100, 10^ and 10^ TeV. The interaction model 
used is the 'reference model'. The curves are analytic calculations, the histograms 
the results of Montecarlo runs. 

Fig. 7a. Probability distribution P(iV^; EQ,E^i^) that a primary proton of energy 
£"0 produces muons above threshold energy Pmin- The points are the results 
of Montecarlo calculation with 6 = 0°, -Emin = 1 TeV, Eq = 10^ TeV. The dashed 
curve is a poissonian distribution with the same average, the solid line is a negative- 
binomial distribution of same average and dispersion, as the montecarlo result. 

Fig. 7b. As in figure 7a, Eq = 10^ TeV. 
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Fig. 8. Function G(x) for two models that yield the same inclusive muon distribu- 
tion. The proton interactions are described by {Kp, P*} = {0.475,0} in model-1 
and {|, 1} in model-2. The other parameters are chosen as in the 'reference model'. 

Fig. 9. Inclusive muon flux calculated with analytic and montecarlo methods 
assuming the primary proton flux 4>o{E) ~ 1.85 E~'^-'^ steepening to oc for 
£; > 3 X 10^ GeV. (see text). 

Fig. 10a. Muon multiphcity distribution 0^(iV^) (^min = 1 TeV and 6^ = 0° ) 
calculated with a montecarlo method assuming a primary cosmic ray flux of protons 
with energy spectrum: (t)o{E) — 1.85 E~'^-'^ steepening to oc E~^ for £■ > 3 x 
10^ GeV . The two set of points refer to two different proton interaction models : 
{Kp, P*} = {0.475, 0} and {Kp, P*} = {§, 1} that result in the same inclusive muon 
distribution. 

Fig. 10b. Ratio of the fluxes 0^(Ar^) of figure 9a, calculated with two different 
models for proton iteractions. 
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